pacbio_ad_pilot_differential_gene_expression_kallisto_counts

1 Method

Gene-level counts were generated from transcript-level isoform counts produced by lr-kallisto (v0.51.1) using the tximport package (v1.32.0) in R. A reference GTF annotation file (Gencode v48) was parsed to create a transcript-to-gene mapping, which enabled aggregation of transcript counts into gene counts.

The gene count matrix was filtered to retain genes with at least 10 counts in at least two samples to ensure robust downstream statistical analysis.

Differential gene expression analysis was conducted using DESeq2 (v1.44.0). Sample metadata included condition (Alzheimer’s disease vs control), and sex as a covariate in the design formula.

Significant genes were defined by an adjusted p-value < 0.05 and a fold-change threshold of |log2 fold-change| ≥ |log2(1.5)|. Volcano plots highlighting significant genes were generated with ggplot2 (v3.5.2).

All analyses were conducted in R (v4.4.0) and tidyverse (v2.0.0) ecosystem packages.

---

# libraries
library(DESeq2)

## Loading required package: S4Vectors

## Loading required package: stats4

## Loading required package: BiocGenerics

## 
## Attaching package: 'BiocGenerics'

## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs

## The following objects are masked from 'package:base':
## 
##     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
##     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
##     get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
##     match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
##     Position, rank, rbind, Reduce, rownames, sapply, setdiff, table,
##     tapply, union, unique, unsplit, which.max, which.min

## 
## Attaching package: 'S4Vectors'

## The following object is masked from 'package:utils':
## 
##     findMatches

## The following objects are masked from 'package:base':
## 
##     expand.grid, I, unname

## Loading required package: IRanges

## Loading required package: GenomicRanges

## Loading required package: GenomeInfoDb

## Loading required package: SummarizedExperiment

## Loading required package: MatrixGenerics

## Loading required package: matrixStats

## 
## Attaching package: 'MatrixGenerics'

## The following objects are masked from 'package:matrixStats':
## 
##     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
##     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
##     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
##     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
##     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
##     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
##     colWeightedMeans, colWeightedMedians, colWeightedSds,
##     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
##     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
##     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
##     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
##     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
##     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
##     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
##     rowWeightedSds, rowWeightedVars

## Loading required package: Biobase

## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.

## 
## Attaching package: 'Biobase'

## The following object is masked from 'package:MatrixGenerics':
## 
##     rowMedians

## The following objects are masked from 'package:matrixStats':
## 
##     anyMissing, rowMedians

library(ggrepel)

## Loading required package: ggplot2

library(rtracklayer)
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.3     ✔ tibble    3.2.1
## ✔ purrr     1.1.0     ✔ tidyr     1.3.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ lubridate::%within%() masks IRanges::%within%()
## ✖ dplyr::collapse()     masks IRanges::collapse()
## ✖ dplyr::combine()      masks Biobase::combine(), BiocGenerics::combine()
## ✖ dplyr::count()        masks matrixStats::count()
## ✖ dplyr::desc()         masks IRanges::desc()
## ✖ tidyr::expand()       masks S4Vectors::expand()
## ✖ dplyr::filter()       masks stats::filter()
## ✖ dplyr::first()        masks S4Vectors::first()
## ✖ dplyr::lag()          masks stats::lag()
## ✖ ggplot2::Position()   masks BiocGenerics::Position(), base::Position()
## ✖ purrr::reduce()       masks GenomicRanges::reduce(), IRanges::reduce()
## ✖ dplyr::rename()       masks S4Vectors::rename()
## ✖ lubridate::second()   masks S4Vectors::second()
## ✖ lubridate::second<-() masks S4Vectors::second<-()
## ✖ dplyr::slice()        masks IRanges::slice()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(tximport)

knitr::opts_chunk$set(echo=TRUE, warning=FALSE, message=FALSE, fig.width=12, fig.height=8, fig.align = "center")

sessionInfo()

## R version 4.4.0 (2024-04-24)
## Platform: x86_64-pc-linux-gnu
## Running under: AlmaLinux 8.10 (Cerulean Leopard)
## 
## Matrix products: default
## BLAS/LAPACK: FlexiBLAS NETLIB;  LAPACK version 3.11.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: America/New_York
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] tximport_1.32.0             lubridate_1.9.3            
##  [3] forcats_1.0.0               stringr_1.5.1              
##  [5] dplyr_1.1.4                 purrr_1.1.0                
##  [7] readr_2.1.5                 tidyr_1.3.1                
##  [9] tibble_3.2.1                tidyverse_2.0.0            
## [11] rtracklayer_1.64.0          ggrepel_0.9.6              
## [13] ggplot2_3.5.2               DESeq2_1.44.0              
## [15] SummarizedExperiment_1.34.0 Biobase_2.64.0             
## [17] MatrixGenerics_1.16.0       matrixStats_1.5.0          
## [19] GenomicRanges_1.56.0        GenomeInfoDb_1.40.0        
## [21] IRanges_2.38.0              S4Vectors_0.42.1           
## [23] BiocGenerics_0.50.0        
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.2.1         farver_2.1.2             Biostrings_2.72.0       
##  [4] bitops_1.0-9             fastmap_1.2.0            RCurl_1.98-1.14         
##  [7] GenomicAlignments_1.40.0 XML_3.99-0.16.1          digest_0.6.37           
## [10] timechange_0.3.0         lifecycle_1.0.4          magrittr_2.0.3          
## [13] compiler_4.4.0           rlang_1.1.6              sass_0.4.10             
## [16] tools_4.4.0              yaml_2.3.10              knitr_1.50              
## [19] S4Arrays_1.4.0           curl_6.4.0               DelayedArray_0.30.1     
## [22] RColorBrewer_1.1-3       abind_1.4-8              BiocParallel_1.38.0     
## [25] withr_3.0.2              grid_4.4.0               colorspace_2.1-1        
## [28] scales_1.4.0             dichromat_2.0-0.1        cli_3.6.5               
## [31] rmarkdown_2.29           crayon_1.5.3             generics_0.1.4          
## [34] rstudioapi_0.16.0        httr_1.4.7               tzdb_0.4.0              
## [37] rjson_0.2.21             cachem_1.1.0             zlibbioc_1.50.0         
## [40] parallel_4.4.0           XVector_0.44.0           restfulr_0.0.15         
## [43] vctrs_0.6.5              Matrix_1.7-0             jsonlite_2.0.0          
## [46] hms_1.1.3                locfit_1.5-9.9           jquerylib_0.1.4         
## [49] glue_1.8.0               codetools_0.2-20         stringi_1.8.3           
## [52] gtable_0.3.6             BiocIO_1.14.0            UCSC.utils_1.0.0        
## [55] pillar_1.11.0            htmltools_0.5.8.1        GenomeInfoDbData_1.2.12 
## [58] R6_2.6.1                 evaluate_1.0.4           lattice_0.22-6          
## [61] Rsamtools_2.20.0         bslib_0.9.0              Rcpp_1.1.0              
## [64] SparseArray_1.4.0        xfun_0.52                pkgconfig_2.0.3

# make t2g file
gtf_file <- readGFF("../../refs/gencode.v48.annotation.gtf")

head(gtf_file)

t2g <- gtf_file %>% dplyr::select(transcript_id, gene_id) %>% 
  na.omit() %>% 
  arrange(gene_id) %>% 
  distinct() %>% 
  as.data.frame()

2 Generate Gene-Level Counts

Using lr-kallisto Isoform Counts

# tximport kallisto transcript counts
kallisto_count_dir <- "../kallisto_counts"
sub_dirs <- list.dirs(kallisto_count_dir, recursive=TRUE)[-1]

files <- file.path(sub_dirs, "abundance.tsv")
names(files) <- stringr::str_sub(basename(sub_dirs),-4)

txi_object <- tximport(files, type = "kallisto", tx2gene = t2g, geneIdCol = gene_id, txIdCol = transcript_id)

gene_counts <- txi_object$counts %>% as.data.frame() %>%
  tibble::rownames_to_column("gene_id")

gene_counts

# read in raw counts matrix
new_col_names <- c("gene_id", "bc01_CTRL", "bc02_CTRL", "bc03_CTRL", "bc04_CTRL", "bc09_AD", "bc10_AD", "bc11_AD", "bc12_AD")

names(gene_counts) <- new_col_names

nrow(gene_counts)

## [1] 78686

3 Filter Counts

# filter counts to keep genes with at least 2 samples with counts >= 10

# all columns to integers
gene_counts[ , -1] <- lapply(gene_counts[ , -1], as.integer)

counts_filtered <- gene_counts[rowSums(gene_counts[-1] >= 10)>= 2, ]

nrow(counts_filtered)

## [1] 18890

counts_filtered

4 Run DESeq2

# create deseq2 object

# counts matrix
counts_matrix <- as.matrix(counts_filtered[-1])
row.names(counts_matrix) <- counts_filtered$gene_id

4.1 Covariate Data for DESeq2

# col data
col_data <- tibble(
  sample_name = colnames(counts_filtered[-1])) %>%
  separate(sample_name,c("sample","condition"),sep="_",remove=FALSE) %>%
  mutate(condition=as.factor(condition))

col_data$sex <- factor(c('male','male','male','female','female','female','male','female'))
#col_data$race <- factor(c('Black','Black', 'White', 'White', 'White','White','Black','Black'))

col_data

4.2 Design Formula

# design matrix and dds object
design <- formula(~ condition + sex)
dds <- DESeqDataSetFromMatrix(
  countData=counts_matrix,
  colData=col_data,
  design=design
)

# run deseq2
dds <- DESeq(dds)
res <- results(dds)

# join gene names

res_df <- res %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene_id") %>%
  arrange(padj) %>%
  inner_join(gtf_file %>% dplyr::select(gene_id, gene_name) %>% distinct(), by = "gene_id")

write.csv(res_df, "kallisto_DGE.csv", row.names = FALSE)

5 Volcano Plot of Differentially Expressed Genes

padj_threshold <- 0.05
log2fc_threshold <- abs(log2(1.5))

# Get genes to label (significant & above FC threshold)
genes_to_label <- res_df %>%
  dplyr::filter(padj <= padj_threshold & abs(log2FoldChange) >= log2fc_threshold) %>%
  pull(gene_name)

# Create gene_label column
res_df <- res_df %>%
  mutate(gene_label = if_else(gene_name %in% genes_to_label, gene_name, NA_character_))

# Make volcano plot
volcano_plot <- res_df %>%
  ggplot(aes(x = log2FoldChange, y = -log10(padj))) +
  geom_point(color = "grey") +
  geom_point(data = subset(res_df, padj <= padj_threshold & log2FoldChange >= log2fc_threshold), color = "red") +
  geom_point(data = subset(res_df, padj <= padj_threshold & log2FoldChange <= -log2fc_threshold), color = "blue") +
  geom_text_repel(
    data = subset(res_df, !is.na(gene_label)),
    aes(label = gene_label),
    max.overlaps = 25,
    box.padding = 0.7,
    point.padding = 0.25
  ) +
  geom_hline(yintercept = -log10(padj_threshold), linetype = "longdash", colour = "grey") +
  geom_vline(xintercept = log2fc_threshold, linetype = "longdash", colour = "red") +
  geom_vline(xintercept = -log2fc_threshold, linetype = "longdash", colour = "blue") +
  labs(
    title = "Volcano Plot of Differentially Expressed Genes",
    x = "Log2 Fold Change",
    y = "-Log10(padj)"
  ) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))

volcano_plot

res_df